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We are developing 3 dimensional simulation codes for coalescing binary neutron stars. A 
code using the maximal slicing condition is obtained. To evaluate the gravitational radiation, 
we implemented a gauge-invariant wave extraction and compared the wave forms with the 
■ metric tensors at the wave zone. The energy spectrum of the waves was also evaluated to 

' investigate the possibility that the excitation of the quasi-normal modes of the black hole, 

which may be formed after the merger of two stars, can be observed. 

y—i ■ 
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1 _— I ' §1. Introduction 

oo : 

We are constructing computer codes on 3D numerical relativityf^^ We report 
O here on the present status of our computer code development to study the fully 

general relativistic evolution of spacetimes and matter. Our main target is to study 
the evolution of coalescing binary neutron stars and the radiation of gravitational 
waves from the merger. Coalescing neutron stars or black hole binaries are the 
^"'i most promising sources of strong gravitational waves for interferometric detectors 

such as TAMA300, LIGO, VIRGO, and GEO600 31 The first detection by these 
detectors may be the waves from a black hole binary because of the larger mass 
of the system. However the astrophysical importance of a coalescing neutron star 
binary is not smaller than a black hole binary. For example, in a certain model 
of a gamma ray bursts j2H3 the central engine is a coalescing neutron star binary. 
Furthermore, a detailed analysis of the detected waves near the merger of two stars 
will give information on the size of a neutron star and then on the equation of state 
of high density matter. 

Gravitational waves from coalescing binaries consist of three phases; (1) the 
inspiral phase, (2) the merging phase and (3) the ringing-down phase. At the inspiral 
phase, the separation between two stars is so large that they may be considered as 
point masses and the luminosity of gravitational waves is so small that the orbit of 
each star may be quasi-stationary. The waves from this phase is called the chirp 
signal and the wave form can be estimated in a post-Newtonian approximation. The 
coalescing binary neutron stars at the merging and the ringing-down phase cannot 
be considered as point masses and the general relativistic effects become large. Thus 
general relativistic numerical simulations are necessary to investigate the evolution 
of the binary at these phases and to predict the radiation of gravitational waves from 
the merger. 
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In the previous code we used the conformal slicing condition in which the metric 
becomes the Schwarzschild one in the outer vacuum region so that the wave extrac- 
tion is easy, while for the long time integration the slicing is not stableP Then we 
started to construct a new code using the maximal slicing condition. Disadvantages 
of the maximal slicing compared with the conformal slicing are 1) the extraction of 
gravitational waves is not straightforward, 2) more CPU hours are required since an 
elliptic partial differential equation should be solved for the lapse function a. For 
the former, Abrahams et-alP* gave the gauge-invariant wave extraction technique 
for axially symmetric, even-parity perturbations in the Schwarzschild metric. We 
extend their technique to non-axially symmetric perturbations both for even- and 
odd-parities. As a result, we found that the gauge dependent modes in a spatial 
metric perturbation is so small at the vacuum exterior region that a good estimate 
of gravitational waves is given easily. Shibata and UryrP have reported gravitational 
waves from the merger of binary neutron stars using a different code with different 
coordinate conditions so that the comparison of our results with theirs is important. 

In the followings, we show recent results of our 3D general relativistic simu- 
lations of coalescence of binary neutron stars and the estimation of gravitational 
waves. In §2 and §3, we describe basic equations for our general relativistic code for 
coalescing binary neutron stars and our coordinate conditions. Numerical methods 
are described in §4 and we show how to extract gauge-invariant gravitational wave 
form in §5. In §6, numerical results of simulations of the coalescing neutron star 
binary are presented and §7 is devoted to the discussions. 



§2. Basic Equations 



In order to follow the time evolution of the space-time and the matter, the 
Einstein equation should be reduced to a system of evolution equations. We use the 
(3+l)-formalism of the Einstein equation. In the (3+l)-formalism of the Einstein 
equation, the line element is given by 



ds 2 



-a 2 dt 2 + 7y (tfe* + P i dt){dx j + (3 j dt), 



(2-1) 



where a, (5 l and 7y are the lapse function, the shift vector and the intrinsic metric 
of 3-space, respectively. The Einstein equation is decomposed to the four constraint 
equations and 12 evolution equations. The Hamiltonian and the momentum con- 
straint equations are 



R + K z 



16vrp H , 



and 



DjiKU-^iK) =8nJi, 
respectively. The evolution equations are 

dtlij = -2aKij + DiPj + Djfli 

and 



(2-2) 
(2-3) 

(2-4) 



d t Kij 



a 



DiDjOt 
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+a (KKij - 2K U K 1 ^ + K U D^ 1 + K Xj Drf 1 + fiDiKy , (2-5) 

where Rij , R = ^ Rij, and K = j^Kij are the 3-D Ricci tensor, its trace, the 
extrinsic curvature tensor, its trace, respectively and Di is the covariant derivative 
associated with 7^ . The quantities p H , J{ and Sij are the energy density, the mo- 
mentum density and the stress tensor, respectively, measured by the observer moving 
along the line normal to the spacelike hypersurface of t = constant. 

In order to give initial data, the constraint equations are solved for given p H and 
Ji. Here we assume that K = and 7^ is conformally flat at t = as 

7i i = 4 7i i , (2-6) 

where 7^' is the flat space metric. Defining the conformal transformation as 

Ky = <t?Kij, Kj = <l>*Kj, K^ = cj) 10 K», p B = 4fip H , Ji = <j> 6 Ji, (2-7) 

we can reduce Ea. l|2-3|) to 

DjK'i = 8irJi, (2-8) 

where Di is the covariant derivative associated with 7«. The traceless extrinsic 
curvature can be expressed as the sum of the transverse traceless part K^ T and the 
longitudinal traceless part as (LW)ij] 



where 



Kij = K$ T + (LW) ij} (2- 



(LW)ij = DiWj + DjWi - \ j {j D l Wi. (2-10) 



Assuming Kj^- = 0, we have the momentum constraint equations fEa. (|2-8|) ) as 



AWi + ^DiDiWj = 8irJi, (2-11) 

where A = &Di. 

The Hamiltonian constraint equation (Eg. (12-2(1 ) becomes 

A<f> = -2n<fr 1 PB-l<fr 7 KijK ij - (2-12) 

o 

To solve the evolution of the metric tensor, we define the following variables as 

= (det( 7ij ))^, (2-13) 

lij = rSij , (2-14) 



(2-15) 
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where 
and 



Kij = 0- 4 (A' ii ) STF , (2-16) 
K = J* ICq (2-17) 
7 ife 7fci = ^-, (2-18) 

(i^-) STF = i + Kji - \iijl kl K k }j . (2-19) 

The indices of K^ are lowered and raised by 7^ and 7^; 

K*j = j ik K kj , & = i jk K\. (2-20) 

We treat 4>, 7y, i 7 *, -fQj and -fT as independent variables in numerical calculations, 
although they are not really independent. In this framework, Eq. ([2-15j) and the 
following equations can be treated as additional constraints, 

det(7y) = 1, (2-21) 

l ik K kl = 0. (2-22) 

Recently some kinds of the reformulation of the Einstein equations in numerical 
relativity have been proposed to obtain numerically stable codes.® Our formulation 
is the simplest one which initiated such researches of reformulations. The motivation 
to use the formulation is that we encountered numerical instability in development 
of a 3-dimensional, fully relativistic numerical simulation code and suppose that nu- 
merical errors in the second derivative of the metric tensor needed to calculate the 
Ricci tensor are likely to cause large errors and the instability. Then we decided to 
compute F l as independent variables and calculate the Ricci tensor using F l (see 
Eas. (|4-2U|) — l|4-24j) ) . Then we found that such modification makes the code dramat- 
ically stable. This formulation is now often cited as BSSN formulation ,®ED but it 
was first introduced by Nakamura, Oohara and Kojimap^ and is continuously used 
in our code development p)ffl£3JI2) 

From Eas. l|2-4|1 and l|2-5|) . we have the evolution equations for these variables: 

dt<f> = ~ {oiK - A/?') , (2-23) 



ft 



tlij 



aK %3 - 0- 4 (A/3i) STF =Aij, (2-24) 



dtKii = <T 4 {a [(R i3 ) STF - 8tt(^) stf ] - (A^a) STF } 



+ a [KKij - 2K u K l j) (2-25) 
+ KuDj/3 1 + k 3l Dif3 l + A U^Kj) - ^KijDtf 1 
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and 

dfK = a 
The quantity F % obeys 



K lj k^ + ^K 2 +47r(p H + S i i ) 



£>*Aa + ^ DiK. (2-26) 
flti* = -A*^-, (2-27) 



since 

da 13 = -i ik i jl da kl = -i ik j jl A kl = -A ij . (2-28) 
We assume that the matter is the perfect fluid, whose stress-energy is given by 

= + pe + p)u^u u + pg^, (2-29) 

where p, e and p are the proper mass density, the specific internal energy and the 
pressure, respectively, and is the four-velocity of the fluid. The energy density 
p H , the momentum density Jj and the stress tensor Sij of the matter measured by 
the normal line observer are, respectively, given by 

Ph = T^iivi Ji = hi^n Tfj,ui Sij = h^hj T^, (2-30) 

where n M is the unit timelike four-vector normal to the spacelike hypersurface and 
is the projection tensor into the hypersurface defined by 

K v = g^u + n^n v . (2-31) 

The relativistic hydrodynamics equations are obtained from the conservation of 
baryon number, V^pu^) = 0, and the energy-momentum conservation law, V U T^ V = 
0. In order to obtain equations similar to the Newtonian hydrodynamics equations, 
we define p N and uf by 

Pn = Viau°p, uf = (2-32) 

respectively, where 7 = det(7jj). Then the equation for the conservation of baryon 
number takes the form 

d t p N +d l [p N V l )=Q, (2-33) 

where 

Vl = i = ^r-P l - ^ 



U u P + p H 
The equation for the momentum conservation is 



dt{p N uf ) + di (p N ufV l ^j = -^Fjadip - ^j(p + p H )<9, ; a 

k T l 



PyaJ J 

+%T- sdilkl + VlJAP 1 . (2-35) 

2(p + Ph) 



The equation for the internal energy conservation becomes 

d t {p N e) + di (p N eV l ) = -pd v {^jau u ) . (2-36) 
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To complete the hydrodynamics equations, we need an equation of state, 

p = p(e,p). (2-37) 

The right-hand side of Ea, (|2-36|) includes the time derivative. For a polytropic 
equation of state, p = (T — l)pe, however, the equation reduces to 

d t (p N e N ) + ch (p N e N V l ) = -p N diV l , (2-38) 

where 

e N = (^jau°) r ~ l e, (2-39) 

p N = (r - l)p N e N = {^au°) r p- (2-40) 

§3. Coordinate conditions 

The choice of the shift vector j3 l and the lapse function a is important because 
the stability of the code largely depends on them and because it is intimately related 
to the extraction of physically relevant information, including gravitational radiation, 
in numerical relativity. Since the right-hand side of Eq. ()2-24|) . Aij is trace- free, the 
determinant of jij is preserved in time and the condition 

DjA ij = (3-1) 

produces the minimal distortion shift vector. This is a good choice of the spatial 
coordinate, but it is too complicated to be solved numerically. Instead we replace 
the covariant derivative in Eq. (3.1) by the partial derivative as 

djlv = 0. (3-2) 

In this condition, we can simply set 

F i = 0, (3-3) 

since Ea. (17771) becomes dtF 1 = and F^t = 0) = 0. 

From Eas. tfPM . (tTTfl) and (J23J), Ea. !^ is reduced to 



0, 



& = -2dj (aK ij ) + ^djdk? + h ij djd k (3 k = 0, (3-4) 



which yields an elliptic equation for the shift vector j3 l : 

V 2 P 1 + U-AP 1 = 2d 3 (aK^ - h? k djd k {3 1 - hi^d 3 d k f3\ (3-5) 
where 

hij = % - Sij and h ij = j ij - S ij (3-6) 

We call this condition as the pseudo-minimal distortion condition. 

As for the slicing condition, we choose the maximal slicing, K = 0. From 
Ea. (|2-26|) . we have an elliptic equation for the lapse function a: 

D m D m a = a (k^K^ + 4vr (p H + S)) . (3-7) 



Numerical Simulation on Coalescing Binary Neutron Stars 



7 



§4. Numerical Methods 

As for the spatial coordinates in numerical calculations, we use a Cartesian 
coordinate system (x,y,z). To obtain initial data, we should solve elliptic partial 
differential equations of Eas. (|2-ll|) and (|2-12|) . The coupled elliptic equations (|2T1I) 
can be reduced to four decoupled Poisson equations: 

A X = GndJi, (4-1) 

AWi = SvrJi - ~d iX , (4-2) 

where \ = ®iWi- These Poisson equations are solved using a pre-conditioned conju- 
gate gradient method™ The boundary conditions for x> Wj and 4> are, respectively, 
given by 

X 2r 3 2r 3 2r 5 \r 4 / 

_ P fc xV 7Pj (7Mjj - Mjj - M^) 3M jk xix k x l ( 1_\ 



and 



M d k x k ^ ( 1 . 



for large r, where 

P,= I Jid 3 x, Mij = I J t x j d 3 x, (4-6) 



M = JSd i x, al k = J Sx k d 6 x (4-7) 

and S 1 is the right-hand side of Ea. (|2T2|) . Since the source term S depends on </>, we 
will find a self-consistent solution of Ea. (|2T2() iteratively. The non-linear iteration 
will usually converge within 10 rounds. 

At each time step, we should also solve the elliptic equation (|3-5|) for the shift 
vector f3 % . The same procedure as for Wi is followed for it. Since the source term 
depends on a self-consistent solution should be found iteratively. The iteration 
will usually converge soon, but it sometimes becomes more than 30 rounds at the 
final stage of a simulation. 

Hydrodynamics equations (|2-33l) . (|2-35|) and (|2-38|) are solved using van Leer's 
schemed The evolution equations of metric Eas. ^FI^i . (J2HU), (|2~23)) and (|2~2H 
are, respectively written as 

d t ^-(3 l d l ^ = -^(aK-d l (3^, (4-8) 



STF* 



(4-9) 
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dfKij — f3 l diKij = 0~ 4 {a [(R, 



\STF 



and 



d t K - (3 l diK = a 
These equations are written as 



„, -8^(^) STF ]-(A^a) STF } 
+ «{KKi, - 2K;, h" , ) + K ll( ), r + Kjtihfi - \K i3 dtf l (4-10) 



K ij K i i + -K 2 +4ir(p H + S i i ) 



d t Q + v%Q = F, 



(4-11) 



(4-12) 



with v = — /3 , which is solved using the CIP (Cubic-Interpolated Pseudoparti- 
cle/Propagation) method. 18 - In the CIP method, both Eq. ()4-12(l and its spatial 
derivatives 

dt(d a Q) + v e d t (d a Q) = -(d l v t ){d e Q) + 8 a F (4-13) 

are solved. This method reduces numerical diffusion during the propagation of Q, 
since the time evolution of Q and its derivatives are traced. Even if the gradient of 
Q becomes very large near the surface of the neutron star, the numerical diffusion is 
small so that errors and the instability does not appear. 

From the Hamiltonian constraint, the conformal factor <j> also obeys 

_ ^5 



A<f> 



16ir PH + KijK 13 



(4-14) 



To obtain <f>, we first make (ft evolve using Ea. (|4-8|) and then calculate the right-hand 
side of Eq. (|4-14|) using new (ft and finally solve Eq. (|4-14|) . 

For the evolution of both matter and metric, we use a two-step algorithm to 
achieve second-order accuracy. That is, to solve the evolution equation such as 



dtQ = F, 



first evolve half a time step hAt, 



Q(t + ±At) = Q(t) + ^At-F(t). 



(4-15) 



(4-16) 



Then F is calculated from quantities at t + \At and finally Q(t + At) is calculated 
from Q(t) and F(t + -At) as 



Q(t + At) = Q(t) + At ■ F(t + -At). 



(4-17) 



We must take special care to calculate the Ricci tensor appearing on the right- 
hand side of Eq. (j4- 1(J|) . With the conformal transformation of Eq. ()2-14|) , the Ricci 
tensor Rij associated with jij is given by 



Rij + Rj<, 



(4-18) 



where 



R% = -20- 1 (DjDrf + jijAcft) + 2(ft~ 2 h{D i( ft)(D j( ft) - %(D k (ft)(D k ^)} (4-19) 
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and Rij is the Ricci tensor associated with 7^. The tensor is given by 



Ra 



l r_ 



7, l k \lli,jk + 7lj,ik - lij,kl) + 1 ki ,l(lki,j + 7kj,i - 1ij,k) 



kl 



2 L 



f*fj fc , (4-20) 



where F % - k is the Christoffel symbol associated with 7^ . The second derivatives of %j 
are replaced by finite differences in numerical calculation but the numerical precision 
of terms such as Ty^j with k^l is not so good, while the degree of precision of 
with k = I is the same as that of the first derivatives. Inaccuracies in 7^ m will cause 
a numerical instability. Since 



7 jk lij,k 



-Hjl jk ,k 



and 
then 



(4-21) 
(4-22) 



R; 



+ 2L 



F k (7ki,j + 7fcj,i _ 7ij,k) 



pA: pi 



On the pseudo- minimal distortion condition, where F l = F{ = 0, we have 



R 



1 



-ki 



-ki 



kl- 



77 17 ,jiu,k-i ,aji,k-j 7%j,ki 



p/c pi 



(4-23) 



(4-24) 



Although the second derivative appears in the term 7 7y jm, it can be written as 

l k %jM = (o~ kl + h kl )h ijM = h ijM + h kl h ijM , (4-25) 

and thus inaccuracies in numerical value of 7^ jy will not be so serious, if both /i^z 
and /li^fci are small. 

Boundary conditions for and 7^ are important, because the boundary is 
not so far from the stars and therefore inappropriate boundary conditions cause 
reflections of the outgoing waves. We apply outgoing conditions for Kij and 7^, in 
our coordinate conditions. The condition can be written as 

H(at - 4> 2 r) 



Q(t,x l 



where r = \J x 2 + y 2 + z 2 . Here H satisfies an advection equation 

d t H + c l diH = 0. 

where 



ax 
r<p 2 



This implies 



a 



d t Q + c^Q 

Ea. (|4-29|) is solved along with Eas.([4-9|) and (|4-ll)|) using the CIP method. 



r<p L 



(4-26) 

(4-27) 
(4-28) 
(4-29) 
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§5. Gauge-Invariant Wave Form Extraction 



In the gauge conditions described in non-wave parts of the perturbations 
decrease as 0(r _1 ) for large r and therefore hij defined by Eq. (|3-6|) cannot be simply 
considered as the gravitational wavesP it includes gauge dependent modes. So that, 
it is necessary to perform a gauge-invariant wave form extraction. Here we apply 
a gauge-invariant wave form extraction technique suggested by Abrahams et.al.P 
which is based on Moncrief's formalism.^ 

Outside of the star, we can consider the total spacetime metric which is 
generated numerically as a sum of a Schwarzschild spacetime and non-spherical per- 
turbation parts: 

a — „i B ) + h(e) i u{o) 

where g^ u is a spherically symmetric metric given by 

/ -iV 2 
A 2 





q (B) 



\ 











R 2 










\ 



(5-1) 



(5-2) 



R 2 sm 2 6 J 



and hiiu and h]^ are even-parity and odd-parity metric perturbations, respectively. 

J 

\ 



They are given by 
(N 2 



and 



Irn 



HllmYlm 

2 ; 



I,( e ) v 
n 0lm 1 lrn.,0 

i + Glm -Qjp 



sym A 2 H 2 i m Yi 



sym 



sym 



sym 



sym 



"'0lm I lm,<t> 
'him 1 lm,<j> 

Yim r Gi m Xi 

m 

^33 



sym 



(5-3) 



E 

lm 



-h 



(o) Yi mt p 
Olm~ 



sm ( 



sym sym 



sym sym 



-h {o) Y 

\ h {o) Xlr 



, sm( 



1 



a sm I 



"Urn 



sm ( 



h 2lm W lrn sin( 



where 



h 



(e) 
33 



r sm 



iV 2 , A 2 , R 2 , Hn m , /ig/m' hii m , Kim, Gkrni K 

and r and Y\ r _ 



sym 
K\ m Y[ m + G; m 

(o) 



-\ h 2im X i™ sin I 



(5-4) 



1 <90 2 



and /i 



W lm 



(5-5) 



(e) 

OZm' "Mm' "<mi ^fcrni '*ojm' "Urn' a,llu ,l, 2Zm 
is the spherical harmonics. The functions Xi m and Wi 

( d 2 8\ 

Xlm = 2 — — - COt 0— Yl m , 



are the functions of t 
are given by 

(5-6) 



Numerical Simulation on Coalescing Binary Neutron Stars 



11 



and 



ill! 



de 2 



„ 9 

cote do 



(J 2 



sin 



Y, 



(5-7) 



respectively. The symbol 'sym' in Eas. ()5-3(l and (|5-4[) indicates the symmetric com- 
ponents. From the linearized theory about perturbations of the Schwarzschild space- 
time, the gauge invariant quantities and are given by 



(o) * 
2lm 



and 



2(/l - 2) ArN 2 k 2 im + Ark 



llm 



A (yl + l-3iV 2 ) 

for the odd and even parity modes, respectively, where 

A = 1(1 + 1), 



and 



k2lm 



K lm + N 2 rG lm>r - 2—h^ m 



H- 



2lm 



1 d 



/ r 



(5-8) 

(5-9) 

(5-10) 
(5-11) 
(5-12) 



2N 2 y/J^dr \N 2 ~ 

The quantities and satisfy the Regge- Wheeler and the Zerilli equations, 
respectively, 

d 2 d 2 



dt 2 dr 2 



(J) 

l in 



0, (/ = e, o) 



(5-13) 



where and are the Regge- Wheeler and the Zerilli potentials.^ The symbol 
r* is the tortoise coordinate defined by r* = r + 2M ln(r/2M — 1). Two independent 
polarizations of gravitational waves h+ and h x are given by 



where 



lYlm 



Lm 



1 



w lr 



(5-14) 



(5-15) 



y/A(A - 2) V sm < 

In numerical calculations, the functions iV 2 (i, r), A 2 (t,r) and R 2 (t,r) of the back- 
ground metric are calculated by performing the following integration over a two- 
sphere of radius r: 



N 2 (t,r) = -±- J g tt df2, 
R\t,r) = ^- [ ( ge9 + -^-) dQ. 



8111 



(5-16) 
(5-17) 
(5-18) 
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where dQ = sin 



The components of the even parity metric perturbations are 
H 2 Ut,r) = grrY^dfl, (5-19) 



Gim(t,r) 



A(A - 2) R 2 



gee 



sin 



2g 9( j > X t 



1 711 



smtf sm( 



K lm (t,r) = l -AG lm + -^J {gee + ^ ) Y lw ) dP 



and 



■'" " f + do. 



KL{t,r) 



A 



suit/ sm( 



(5-20) 
(5-21) 

(5-22) 



where * denotes the complex conjugate. The components of odd parity metric per- 
turbations are 



1 

A 



lm,4> gr4> -rrif \ | 

g r e— - - —^Yim.e I d ^ 



sm I 



sm( 



and 



A(A - 2) 



sm 



X. 



lm 



^ge<t> 



W, 



sm U sm ( 



d(2. 



(5-23) 



(5-24) 



The metric tensor g rr , ggg etc. on the spherical coordinate system appearing 
in Eqs. (|5-17j) - (|5-24|) are calculated from g xx , g yy etc. on the Cartesian coordinate 
system, such as 



g rr = sin 2 (g xx cos 2 + 2g xy sin<f> cos <f> + g yy sin 2 ) 
+gzz cos 2 9 + 2 sin6* cos 9 (g yz sm<fi + g zx cos (f> ) . 

We need angular integrals over spheres for constant r, such as 



(5-25) 



F(ro)= f(x,y,z)df2= f(r ,9,<t>)df2, (5-26) 

J r=ro J 

while we have the physical quantities only at Cartesian grid points. We need in- 
terpolation to obtain the values of f(ro,9,4>) from f(x,y,z) at the grid points. It 
is, however, not easy to fully parallelize the procedure on a parallel computer with 
distributed memory. We therefore rewrite Eq. (|5-2fi|) as a volume integral, namely, 



F(r ) 



f(x, y, z)6(r — ro) d 3 x = lim 



1 



a^o \/TTar, 



'o 



F(x,y,z)e- < - r - ro ^^ 2 d 3 x, (5-27) 



where r = \J x 2 + y 2 + z 2 . Numerical integral with a = Ax/2 gives a good value to 
Ea. (|5-27|) . where Ax is the separation between grid points. It is favorable both in 
the accuracy and the speed for parallel computers. 
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§6. Numerical results 

We performed numerical simulations for coalescing binary neutron stars and 
evaluated the gravitational waves. We used 475 x 475 x 238 Cartesian grid assuming 
the symmetry with respect to the equatorial plane, which require memory of about 
80GBytes. The separation between grid points is Ax = Ay = Az = 1M , i.e., equal 
spacing. Here we use the units of G = c = 1. As the initial condition, we put two 
spherical stars of rest mass 1.5M® and radius 7.7M = 11.6km. The separation 
between the center of each mass is 3O.2M = 46.2km. As for an equation of state, 
we use the 7 = 2 polytropic equation of state. The initial rotational velocity is given 
so that the circulation of the system vanishes approximately as 

V a (r) = fixr-nx{r-r a ), (a = 1,2) (6-1) 

where 17 is the orbital angular velocity and r a is the location of the center of each 
star. Now we set Q = O.O1O/M , then the angular momentum Jo becomes 6.7Mq. 
The total ADM mass M A dm of the system is 2.76M and q = Jo/M ADM = 0.89. 
Figures ^ shows the evolution of the density p N on the x-y plane. The stars 
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Fig. 2. Plots rh+ iX along 2-axis at r = 110 ~ 
140 A/© as a function of t — r. 
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Fig. 3. Plots r 3 /i+, x as a function of t. 



start to coalesce at approximately t = 0.5msec and an almost axisymmetric star is 
formed by t = 1.8msec. 

FigureEJshows the gravitational wave forms rh + and rh x on the z-axis evaluated 
at r = 110, 120, 130 and 14OM as functions of the retarded time t — r. The lines 
of rh+(t — r) and rh x (t — r) estimated at r = 110 ~ 14OM for t — r > coincide 
with each other but they don't for t — r < 0. Figure |3] shows r 3 /i as a function of t. 
It reveals that h includes a non-wave mode proportional to r~ 3 , which corresponds 
to the quadrupole part in the Newtonian potential of the background metric. It 
decreases fast as the merger of stars proceeds. 

Therefore we eliminate the non-wave mode from /i+ and h x using Fourier trans- 
formation as follows: 

• Assuming that total waves are expressed as a sum of wave parts F(t — r) jr and 
non-wave parts G(t)/r 3 , 

\ F(t-r) G(t) , , 

h(t,r)= K r ' +-^. (6-2) 

• Fourier components of h(t, r) are written as 

K{r) = F w (r) + -gG w (r). (6-3) 

F ^^S F( - t)e ~ iujtdt (64) 

and 

G u = ^- I G(t)e-^dt. (6-5) 



(6-6) 



where 



2?r . 

From the values of h u {r) in different radial coordinates r\ and ri-, F u can be 
given by 

u r^e~ 1UJT2 y>2g— iwn 
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Fig. 4. Wave forms rh +tX along z-axis as a func- 
tion of t — r. The curves are averages of rh+ tX 
estimated at r = 110 ~ 2OOM0 and error 
bars denote 2a. 



10 - 




-0.5 0.5 1 

t—r (msec) 



Fig. 5. The comparison h+, x shown in FielH 
(solid lines)and h+ lX defined by Eqs. 16-91 1 
and 16T0I (dashed lines). 



Fourier components h w (r) of gravitational waves are given by 

h w (r k ) = F u , k = 1, 2. (6-7) 

By inverse Fourier transformation of h u (r), we can get the gravitational waves, 
which do not include non-wave modes, 

\U)t 



h+(t, r) - ih x (t, r) = J K{r)e lu)t do;. (6-8) 

The resultant wave form is shown in Fig@] The curves represent the average of /i+ 
and h x calculated at r = 110, 120, • • • 2OOM0 and twice the dispersion 2a is shown 
as error bars. 

Figure [5] compares h + and h x given by Eq.JfFHJ) with the metric perturbation 
h + and h x defined by 

1 
2 
and 

h x = h xy , (6-10) 



h+ = - (h xx - h yy ) , (6-9) 



respectively, where hij is defined by Ea. (|3-6j) , Although h + and h x include gauge 
dependent mode, they almost coincide with the gauge independent wave /i+ and h x . 
It shows that the gauge mode in h + and h x is small, while it may be only accidental. 
However, we note that the pseudo-minimal distortion condition Eq. ()3-2(l guarantee 
Ylj fyhij = 0, and thus h+ and h x are transverse-traceless. Then we can argue "the 
energy density of the gravitational waves" as 
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Fig. 6. The propagation of the gravitational waves. "The energy densities of gravitational waves" 
r 2 p„ w on the x-y plane are shown as gray scale figures. Time in units of milliseconds is shown. 



The propagation of r 2 /)Q\v is shown in Fig|BJ A spiral pattern appears, which can 
be explained naively by the quadrupole wave pattern given by 



(6-12) 
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Fig. 7. The luminosity dEQ-yy/dt (in units of Mqc 2 /sec) and total energy AEqw (in units of 
10 _3 Mqc 2 ) emitted as gravitational waves up to t. The solid and dashed lines plot the values 
calculated using Ea, 16T3> and Ea l6T5l . respectively. They are evaluated at r = 2OOA/0. 



On the x-y plane, where 6 = tt/2, pew is constant along the spiral of r + <fr/f2 = 
constant. 

The gauge invariant luminosity d-Ecw/di and the total energy AEq^j of the 
gravitational waves can be calculated as 

^fr = ^£ (K^f + K'M 2 ) (6-i3) 



2,m 

and 



AE GW = J^dt, (6-14) 

respectively. The luminosity and the total energy can be also estimated from Eq. (|6-ll|) 

as 



dt 

and 



2 

/ 

r— >oo 



r 2 powdi7 (6-15) 



AE GW = J^dt, (6-16) 

which is not gauge invariant. The luminosity and the total energy emitted as gravi- 
tational waves up to t calculated at r = 200M Q using Eas. 1)6-13(1 and (|6T4(I as well 
as Eqs. (|6-15|) and (|6-16|) are plotted as functions of t — r in Fig0 which shows that 
Ea. (|6-ll|) gives a good estimate of the energy density of the gravitational waves as 
expected from Fig|SJ 
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Fig. 8. Lapse function a (solid lines) and density p N (dashed lines) along x-axis at t=0(left) and 
1.72msec (right). 

§7. Conclusion 



We showed a stable code using the pseudo-minimal distortion condition and the 
maximal slicing condition for a coalescing neutron star binary. 

We were able to extract the wave form of the gravitational radiation using the 
gauge-invariant extraction techniques. The amplitude of the gravitational waves is 



ft~lx 1(T 21 ( —?—) . 

As shown in Fig(7J the total energy emitted is 

AE GW (t < 3msec) « 2 x l(T 3 M c 2 «3x 10 51 erg » 0.1% of M tot c 2 . 
The angular momentum lost by the gravitational waves are given by 



AJ 



1 

32tt 



/E(K£^2| + h^^ (0) 

Lra 



dt, 



(7-1) 



(7-2) 



(7-3) 



which is 7 x (GMq/c) ~ 1% of the initial angular momentum of the system. 

We have not searched an apparent horizon to see if a black hole is formed, Instead 
we plot a and p N along the x-axis at t=0 and 1.72msec in FigEJ The gradient of a 
near the surface of the merged star becomes large and the value gets less than 0.3. 
This suggests the formation of a black hole. 

If the black hole is formed after the merger of two neutron stars, quasi-normal 
modes of the black hole are excited. To investigate a possibility that the excitation 
of the quasi- normal modes can be seen by the numerically calculated waves, we 
evaluated the energy spectrum of the gravitational waves, which is given by 



dto 



9.-7T ^ 



32vr 



Lm 



+ 



(7-4) 
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Fig. 9. The energy spectrum of the gravita- 
tional waves. The curves are averages of 
dE/duj estimated at r = 110 ~ 200M Q and 
error bars denote 2a. 



where ^itL,{r) is the Fourier transfor- 
mation of &irr}(t,r), and is shown in 
Fig. |UJ The frequencies of the quasi- 
normal modes depend on the angular 
momentum of the black hole, but the 
fundamental frequency of I = 2 for a 
Schwarzschild black hole of mass 2.8M 
is uj = 25 msec -1 . A peak near this 
frequency appears in Fig. |5J Unfortu- 
nately, however, rotating angular fre- 
quency just when the merger of the stars 
finishes is 12 ~ 15 msec - and thus 
they will radiate the waves of frequency 
near uj = 25 msec -1 . So it is not clear 
whether this peak corresponds to the 
emission of the quasi-normal mode of 
the formed black hole. 

Our results must be compared with 
the model H-l of Shibata and Uryu.^ 

The wave form and the luminosity of the gravitational waves are qualitatively con- 
sistent, while our values of the radiated energy and angular momentum are smaller 
than theirs. It caused by the difference of the initial data; two stars merge more 
quickly in our calculation since they are not in quasi-equilibrium at the initial time. 
Nevertheless, the ratio of (Z\£?gw/-^adm) to (AJ/Jq) agrees with their discussion. 

Finally, we must mention the accuracy of our code. Essential parts of code 
tests are summarized in Oohara and NakamuraP We need not solve Ea. (|2-26|) for 
K since we use the maximal slicing, where K = 0, but it is solved to monitor the 
numerical precision. The value of K is kept as small as 0.1% of a typical value of 
Kij but it becomes large from the numerical boundary region in the final stage of 
simulation. It is due to a small reflection of gravitational waves at the numerical 
boundary. Since we use the pseudo-minimal distortion condition, F l = 7^ must be 
zero. It is satisfied within a few % of the derivative of 7 1 - 7 , while the error grows in 
the same way as K in final stage of simulation. The total ADM mass is conserved 
within 10% error up to t < 1.7msec. When two stars merge and a single compact 
object is formed, the conservation get worse, since our grid is too coarse to represent 
such a compact object. Since it is likely to be a black hole by that time, however, 
the global characteristics of gravitational radiation should be hardly affected. 

We found that the gravitational waves are in small part reflected at the numerical 
boundary and reflected waves grow up gradually after t ~ 2.5msec. It is because 
the boundary is still too close. In reality, if the distance of the numerical boundary 
is halved, the growth of the reflected waves near the boundary will be visible by 
t = 1.4msec in the distribution of the energy density of the gravitational waves P GW) 
while it is not yet visible in FigHO This problem must be overcome if the numerical 
boundary is located farther. In addition, a finer grid is required to investigate the 
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formation of a black hole. To perform simulations with a finer and larger-sized 
grid, we must improve schemes for solving elliptic partial differential equations and 
evolution equations. The most CPU hours are consumed for solving elliptic equations 
for /? (Eqs.flO!) and (EHO, « (Eq.(|S3J) and 4> (Ea. QFTty ). and thus the reduction of 
CPU hours required to solve these equations is effective. The present code requires 
memory of 80GBytes and 100 CPU hours with a 475 x 475 x 238 grid. The size is not 
restricted by the memory but by CPU hours. Therefore, larger-scale simulations can 
be performed according to the speedup of the code. We are still improving the code, 
including implementation of new schemes and an apparent horizon finder as well 
as reformulation of the Einstein equations, some results of which will be presented 
soon. 
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